a=[1,1,1,0,-1,0,0,0,0,0,0,0,0;
  -1,0,0,1,0,0,0,0,-1,0,0,0,0;
  0,-1,0,0,1,1,0,0,0,0,0,0,0;
  0,0,-1,0,0,0,1,1,0,0,-1,0,0;
  0,0,0,-1,0,0,-1,0,1,1,0,0,-1;
  0,0,0,0,0,-1,0,-1,0,0,1,1,0;
  0,0,0,0,0,0,0,0,0,-1,0,-1,1];
b1=[4;-4;0;0;0;0;0];
b2=[4;0;-4;0;0;0;0];
b3=[0;-4;4;0;0;0;0];
b4=[4;0;0;0;0;0;-4];
c=[zeros(52,1);ones(13,1)];
F=[eye(13),eye(13),eye(13),eye(13)];
Als=[1*F,-1*eye(13);
     3*F,-1*eye(13);
    10*F,-1*eye(13);
    70*F,-1*eye(13);
   500*F,-1*eye(13);
  5000*F,-1*eye(13);
   -1*eye(13),zeros(13),zeros(13),zeros(13),zeros(13,13);
    zeros(13),-1*eye(13),zeros(13),zeros(13),zeros(13,13)
    zeros(13),zeros(13),-1*eye(13),zeros(13),zeros(13,13)
    zeros(13),zeros(13),zeros(13),-1*eye(13),zeros(13,13)];
cap=5*ones(13,1);
bls=[zeros(13,1);2/3*cap;16/3*cap;178/3*cap;1468/3*cap;16318/3*cap;zeros(52,1)];
Aeq=[a,zeros(7,13),zeros(7,13),zeros(7,13),zeros(7,13);
    zeros(7,13),a,zeros(7,13),zeros(7,13),zeros(7,13);
   zeros(7,13),zeros(7,13),a,zeros(7,13),zeros(7,13);
    zeros(7,13),zeros(7,13),zeros(7,13),a,zeros(7,13)];
beq=[b1;b2;b3;b4];
[x,fval]=linprog(c,Als,bls,Aeq,beq,[],[]);
reslut=[x(1:13,1:1),x(14:26,1:1),x(27:39,1:1),x(40:52,1:1)]
opt=fval